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Abstract 


Elastic  scattering  by  structures  embedded  in  an  ocean  acoustic  waveguide  is  a  problem  of 
paramount  importance.  Apart  from  methods  such  as  the  T-matrix  method,  numerical  methods  such 
as  the  finite  element  method  and  virtual  source  approach  have  been  developed  to  study  this 
problem.  The  present  study  is  concerned  with  the  development  of  a  finite  element  model  for  elastic 
scattering  by  objects  in  depth  and  range  dependent  ocean  waveguides.  The  solution  to  the  three 
dimensional  Helmholtz  equation  in  a  cylindrical  coordinate  system  has  been  expanded  in  a  Fourier 
series  in  the  azimuthal  coordinate.  The  resulting  two  dimensional  Helmholtz  equation  for  each 
Fourier  harmonic  has  been  solved  using  a  finite  element  model.  The  acoustic  pressure  field  in  the 
waveguide  due  to  a  point  source  insonifying  an  elastic  object  in  the  waveguide  has  been 
decomposed  in  to  rigid  scattering  and  radiation  components.  An  assumed  mode  model  has  been 
chosen  to  formulate  the  radiated  field  due  to  structural  vibration.  The  radiation  loading  has  been 
identified  in  terms  of  added  mass  and  radiation  damping  and  the  associated  waveguide  models 
have  been  discussed.  The  FE  model  for  the  waveguide  requires  the  depth  eigensolution  at  the 
truncation  boundary  in  order  to  impose  absorbing  boundary  conditions.  A  Rayleigh-Ritz 
approximation  has  been  developed  to  obtain  the  depth  eigensolution.  The  depth  modes  have  a 
compact  fonn  consisting  of  a  linear  combination  of  isovelocity  modes.  The  research  work  is  being 
pursued  employing  the  FORTAN90  code  developed  implementing  the  proposed  finite  element 
model  for  the  waveguide  scattering  problem. 


DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


Chapter  1 
Introduction 


1.1  General 

Elastic  scattering  by  structures  submerged  in  an  ocean  waveguide  is  a  problem  of  paramount 
importance.  The  waveguide,  by  virtue  of  its  boundaries  and  possible  heterogeneity  of  the  acoustic 
medium,  would  influence  the  elasto-acoustic  behavior  of  the  structure  and  far-field  acoustics  in 
the  waveguide.  A  study  of  this  problem  requires  a  structural  dynamic  formulation,  a  formulation 
dealing  with  acoustic  propagation  in  the  waveguide  and  appropriate  interface  conditions  at  the 
fluid-structure/body  interface.  The  present  study  is  primarily  concerned  with  waveguide  acoustics, 
but  the  motivation  for  the  study  is  established  by  considering  a  coupled  formulation  involving 
structural  dynamic  and  fluid  dynamic  equations,  which  is  familiarly  known  as  the  fluid-structure 
interaction  (FSI)  problem. 

1.2  Literature  Review 

Rigid  body  scattering  and  elastic  scattering  from  objects  in  an  unbounded  domain  have  been 
studied  by  analytical  methods  for  simple  geometries.  Powerful  numerical  methods  such  as  the 
transition  matrix  (T-matrix)  method,  boundary  element  method  (BEM),  and  finite  element  method 
(FEM)  have  been  developed  over  the  years  to  study  the  scattering  problem  involving  complex 
geometries.  In  this  context,  an  important  component  of  domain  discretization  method  such  as  the 
FEM  is  the  use  of  appropriate  nonreflecting  or  absorbing  boundary  condition  at  a  boundary  that 
truncates  an  otherwise  infinite  domain.  Scattering  analysis  of  bodies  in  an  ocean  waveguide 
presents  additional  difficulties  because  of  the  sea  surface  and  the  seafloor,  apart  from  possible 
stratification  of  acoustic  properties  of  the  medium. 

In  most  studies,  the  wave-theoretical  fonnulation  has  been  employed  assuming  the  waveguide  to 
be  cylindrically  symmetric,  which  then  restricts  the  scatterer  geometry  to  be  axisymmetric. 
Hackman  and  Sammelmann  (1986)  extended  the  T-matrix  fonnulation  to  scattering  in 
inhomogeneous  waveguides.  The  analysis  is  done  in  two  steps.  The  first  one  involves  the  analysis 
of  an  empty  waveguide  (i.e.  without  the  scatterer  in  it)  with  a  point  source  for  which  a  normal 
mode  approach  has  been  adopted.  The  second  problem  considers  the  same  waveguide  with  the 
scatterer  located  in  it  for  which  the  T-matrix  method  has  been  employed.  They  have  observed 
multiple  scattering  of  acoustic  waves  by  the  scatterer  and  waveguide  boundaries,  and  the 
inhomogeneous  layers  in  the  waveguide  may  be  important  in  shallow  water  waveguides.  They 
have  further  observed  that  even  in  long  range  propagation,  multiple  scattering  effects  might  be 
important  if  there  is  appreciable  sound-speed  gradient  in  the  vicinity  of  the  scatterer. 

Ingenito  (1987)  presented  a  method  for  scattering  analysis  of  axisymmetric  bodies  in  a  horizontally 
stratified  waveguide.  He  used  the  waveguide  Green’s  functions.  Thus  this  model  fully  accounts 
for  multiple  reflections  at  the  waveguide  boundaries,  but  it  does  not  provide  for  multiple  reflections 
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at  the  scatterer  surface  because  of  the  use  of  free  space  scattering  function.  For  this  reason,  this 
model  is  called  a  single-scatter  model.  Ingenito  (1987)  presented  results  for  a  spherical  scatterer 
and  outlined  the  formulation  for  a  non-spherical  body.  He  has  noted  that  the  method  may  fail  for 
depth-dependent  sound  speed  profdes  because  of  the  presence  of  evanescent  modes. 

Hackman  and  Sammelmann  (1988)  presented  multiple  scattering  analysis  for  a  target  in  a 
waveguide.  Unlike  Ingenito  (1987),  these  authors  have  used  free-space  Green’s  functions.  They 
have  noted  that  rescattering  from  the  object/target  located  in  a  shallow  water  waveguide  could  be 
important  if  there  is  significant  amount  of  radiation  from  it  in  the  vertical  direction. 

Collins  and  Werby  (1988)  presented  a  simple  and  efficient  method  using  the  parabolic  equation 
model  for  scattering  in  a  range-dependent  waveguide. 

Nonnan  and  Werby  (1991)  presented  a  T-matrix  model,  which  implies  a  single  scatter 
approximation  similar  to  Ingenito.  They  have  shown  through  examples  that  the  scattering  object 
is  correctly  coupled  to  the  waveguide.  An  object  is  said  to  be  correctly  coupled  to  the  waveguide 
if  the  scattered  field  displays  the  same  characteristics  as  the  field  produced  by  a  point  source. 

Perkins  et  al.  (1992)  extended  Ingenito’ s  modal  approach  to  weakly  range-dependent  problems 
using  the  adiabatic  approximation. 

Makris  (1998)  presented  a  spectral  fonnulation  for  studying  the  scattering  characteristics  of  a  non¬ 
compact  (i.e.,  1  <  ka,  the  scattering  parameter)  object  in  a  shallow  water  waveguide,  by 
generalizing  the  single-scatter  modal  approach  of  Ingenito. 

The  foregoing  works  employ  one  or  more  of  the  following  restrictions  to  facilitate  the  scattering 
fonnulation:  i)  the  waveguide  is  range-independent  or  at  best  weakly  range  dependent  and  ii) 
single-scattering  which  ignores  multiple  scattering  from  the  object;  iii)  The  object  is  far  removed 
from  the  boundaries  of  the  waveguide  so  that  free  space  scattering  functions  of  the  object  can  be 
used;  and  iv)  the  object  is  sufficiently  far  from  the  source  and  the  receiver.  The  coupled  mode 
method  of  Evans  (1983),  and  full-fledged  numerical  methods  such  as  the  FEM  and  virtual  source 
approach  (VS A)  can  handle  the  general  case  of  scattering  in  range-dependent  waveguides. 
Zampolli  et  al.  (2007)  developed  a  finite  element  model  for  scattering  by  an  axisymmetric  structure 
immersed  in  a  range  and  depth  dependent  cylindrically  symmetric  waveguide.  They  employed 
perfectly  matched  layers  (PML)  at  the  truncation  boundaries  to  simulate  absorbing  boundary 
conditions.  In  the  VS  A,  which  has  been  briefly  reviewed  in  Jensen  et  al  (20 11)  with  key  references, 
several  point  sources  of  unknown  intensity  are  distributed  over  the  scatterer  and  the  source 
intensities  are  obtained  by  solving  a  set  of  global  equations.  The  main  thrust  of  the  present  study 
is  to  formulate  a  3-D  finite  element  model  for  range  and  depth  dependent  cylindrically  symmetric 
waveguides  with  axisymmetric  objects  embedded  in  them. 


DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


1.3  Aim  and  Scope  of  the  Study 


The  primary  aim  of  this  study  is  to  develop  a  finite  element  model  for  elastic  scattering  by 
axisymmetric  bodies  submerged  in  a  heterogeneous  ocean  acoustic  waveguide.  A  3-D  fonnulation 
for  nonaxisymmetric  scattering  and  radiation  in  a  cylindrically  symmetric  waveguide  is  considered 
wherein  the  azimuthal  variation  in  the  acoustic  field  is  represented  using  a  Fourier  series  and  the 
Helmholtz  equation  for  each  of  the  hannonics  is  solved  using  a  finite  element  model.  As  a  part  of 
the  present  work,  a  Rayleigh-Ritz  approximation  for  the  depth  eigenproblem  has  been  developed. 
The  depth  eigensolution  obtained  using  this  model  has  been  employed  at  the  truncation  boundary 
of  the  FE  model  for  the  waveguide. 

For  convenience,  the  elastic  scattering  problem  is  decomposed  into  rigid  scattering  and  radiation 
problems.  In  order  to  bring  out  the  relevance  of  waveguide  radiation  analysis,  a  compact  dynamic 
model  for  a  structure  vibrating  in  contact  with  a  fluid  is  discussed  in  Chapter  2.  3-D  acoustic 
waveguide  propagation  models  corresponding  to  incident  wave,  scattering  and  radiation  problems 
are  discussed  in  Chapter  3.  A  finite  element  model  for  solving  a  typical  waveguide  problem  is 
presented  in  Chapter  4.  A  Rayleigh-Ritz  approximation  is  presented  for  the  solution  of  the  depth 
eigenproblem  in  Chapter  5.  The  eigensolution,  which  is  in  a  compact  form,  has  been  used  in  the 
finite  element  model.  Salient  observations  are  made  in  Chapter  6. 
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Chapter  2 

Acoustic  Radiation  Loading  on  Structures  Immersed  in  an 
Ocean  Acoustic  Waveguide 


2.1  Introduction 

Fluid-structure  interaction  (FSI)  fonnulations  involving  vibrations  of  structures,  such  as  plate  and 
shell  structures,  in  contact  with  a  fluid  are  coupled  in  nature.  A  typical  FSI  formulation  involves 
structural  dynamic  and  fluid  dynamic  equations,  which,  by  virtue  of  coupling  between  them,  have 
mutual  influence.  In  this  context,  the  acoustic  fluid  model  (i.e.,  a  compressible  inviscid  fluid 
model,  Junger  and  Feit;  1972)  is  considered  adequate,  although  viscous  compressible  fluid  models 
are  now  feasible  with  the  emergence  of  powerful  CFD  codes,  but  such  an  option  would  be 
computationally  voluminous.  FSI  problems  are  generally  solved  using  numerical  methods  such 
as  the  BEM  and  FEM.  A  compact  analysis  technique  for  relatively  simple  FSI  problems  consists 
of  the  use  of  an  assumed  mode  expansion  for  the  structural  displacements  within  the  framework 
of  a  variational  formulation,  which  results  in  a  set  of  coupled  discrete  equations  in  the  time  domain. 
Time  harmonic  FSI  problems  have  an  interesting  feature  namely,  the  hydrodynamic/acoustic 
problem  can  be  solved  first  for  each  of  the  assumed  modes  having  unit  amplitude  (Junger  and  Feit, 
1972).  The  results  of  the  acoustic  problem  can  then  be  incorporated  in  the  structural  dynamic 
equations  in  the  form  of  modal  fluid  loads  with  unknown  amplitudes.  The  solution  of  the  structural 
dynamic  equation  at  this  stage  essentially  solves  the  problem  fully,  i.e.,  both  the  structural  and  fluid 
response  quantities  are  determined  at  this  stage.  This  approach  may  be  looked  upon  as  a  sequential 
analysis  procedure  applied  to  an  otherwise  coupled  problem  which  results  in  great  computational 
economy.  Such  an  approach  is  also  widely  employed  in  the  study  of  rigid  body  motions  of  floating 
systems  such  as  ships  subjected  to  small  amplitude  monochromatic  water  waves. 

A  brief  outline  of  the  FSI  fonnulation  is  first  presented  and  then  specialized  to  the  case  of  harmonic 
response,  in  the  following  sections.  The  fluid  loading  on  the  structure  due  to  its  motion  may  be 
obtained  by  solving  the  acoustic  waveguide  propagation  problem,  for  a  set  of  assumed  modes  of 
unit  amplitude.  The  fluid  loading  is  often  referred  to  as  the  acoustic  radiation  loading.  This 
provides  the  motivation  for  the  study  of  acoustic  radiation  from  bodies  located  in  an  ocean 
waveguide. 

2.2  Formulation  of  FSI  Problems 

Consider  small  amplitude,  elastic  vibrations  of  a  shell  in  contact  with  an  acoustic  fluid.  Let 
jwj  =  (u,v,  vt’}7  denote  the  time-dependent  meridional,  circumferential  and  radial  displacements 

of  the  shell.  Let  the  generalized  strain  array  be  denoted  by  {f}  and  the  corresponding  (symmetric) 
elasticity  matrix  by[f9] .  The  strain  array  may  be  related  to  the  displacement  vector  as 
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{*}  =  [!]{£}  (2.1) 

where  the  operator  matrix  [/,] ,  which  denotes  the  generalized  strain-displacement  relation,  may 

readily  be  identified  from  the  shell  theory  adopted  (see  Junger  and  Feit,  1972).  Then,  the  strain 
energy  of  the  shell  is  given  by 

0  =  (2.2) 

^  s 

where  S  denotes  the  surface  of  the  shell.  The  kinetic  energy  may  be  written,  after  neglecting  the 
rotary  inertia  contribution,  as 


(2.3) 


where  m  denotes  the  mass  per  unit  surface  area  of  the  shell  and  a  (.)  indicates  differentiation  with 
respect  to  time  t. 

As  an  approximation,  let  the  displacement  field  be  expressed  as  a  linear  combination  of  a  number 
of  appropriate  assumed  modes,  say  U1,U2,  . . . ,  Un] ,  for  u  ,  VX,V2  ,  •■■,Vn2,  for  v ,  and  WX,W2  ,  . . . , 
Wn 3 ,  for  w  all  being  functions  of  the  spatial  variables  only.  Then,  the  displacement  components 
may  be  written  as 


u 

\uuu2>. 

ir 

o 

p 

V 

>  = 

0,0,... 

.  Vuv2,. 

•••V„2  ,o,o,.... 
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. wvw2,. 

w 

‘“nn3  J 

(avd2>. 

. d„JJ2... 

..,bn2,c1,c2>... 

-cb3) 

(2.4a) 


or 


{if}  =  [A]{J}  (2.4b) 

where d. \t),j  =  1,2, (/?,  +n2  +n3)  denote  the  unknown  generalized  coordinates  which  are 

functions  of  time  t.  For  convenience,  define  an  array  {  N  j  such  that 

{w}r  =  {U„U2,  .M„I,Vl,V2,...V„2,ir],  W2,..W„)  (2.4c) 

which  will  be  used  later.  In  the  view  of  Eq.  (2.4b),  the  strain  array  in  Eq.  (2. 1)  may  be  written  as 
{«}  =  [L][jq{d}  =  [«]{</}  (2.5) 
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In  view  of  Eq.  (2.5),  the  strain  energy  and  kinetic  energy  of  the  shell  respectively  in  Eqs.  (2.2)  and 
2.3  may  be  written  as 

d^dS  (2.6) 

2  s 

T=\\m{^\  (2'7) 

Let  the  excitation/dynamic  load,  including  fluid  dynamic  loading  at  the  fluid-solid  interface,  acting 
on  the  structure  be  denoted  as 


{/}  =  (/,./„./„)'  (2-8) 

where  the  three  components  on  the  r.h.s.  correspond  to  the  three  coordinate  directions.  The  virtual 
work  done  by  the  dynamic  load  may  be  written  as 

SW(t)  =  | [SU}T  \}ys  =  | { Sd}T  [TVf  [  f  }dS  (2.9a) 

or,  using  the  notation  in  Eq.  (2.4c), 


SW  =  Z  J  fuNjSdjdS  +  £  J  ANjSdjdS  +  £  J  fJtjMjdS 

j=l  S  y=«i+lS  j=n,+n2+l  S 

which  may  be  written  symbolically  as 


(2.9b) 


+fl2  +«3 

sw=  £  Qj(t)Sdj  (2.9c) 

7=1 

where  (r)  denotes  the  generalized  load  associated  with  d. .  The  fluid  dynamic  loading,  which  is 
part  of  the  generalized  load,  will  be  discussed  in  Sec.  2.3  dealing  with  hannonic  response. 


The  Lagrange  equations  of  motion  of  the  shell  may  be  written,  using  Eqs.  (2.6),  (2.7)  and  (2.9),  as 


d 

dt 


(  3 

ST 

'  dU 

W-J 

=  Q, 


(2.10) 


Upon  using  the  expressions  for  K.E.  and  P.E.  given  in  Eqs.  (2.6)  and  (2.7),  the  matrix  equation  of 
motion  is  obtained  from  Eq.  (2.10)  as, 

[M]{rf)+[/f]{</}={e}  (2Ji) 

where  the  symmetric  stiffness  matrix  is  given  by 

[Z]  =  j[5f  [Z)][5]rfS  (2.12a) 


and  the  symmetric  mass  matrix  is  given  by 

[M]  =  J/w[JVf 


(2.12b) 
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Note  that  Eqs.  (2.12a)  and  (2.12b)  are  valid  for  time-invariant  heterogeneous  elastic  and  mass 
properties.  The  surface  integrals  above  may  in  general  be  evaluated  using  numerical  quadrature. 

It  may  be  noted  that  the  generalized  fluid  load  component  in  |g|  (see  Eq.  (2.11))  is  unknown, 

which  can  only  be  obtained  by  solving  the  fluid  dynamic  equation.  But  such  a  solution  requires 
constraints  imposed  by  the  fluid  on  structural  motion  at  the  fluid-solid  interface  which  is  also 
unknown.  Thus,  the  structural  dynamic  equations  in  (2.1 1)  and  the  fluid  dynamic  equations  for  the 
waveguide  (discussed  in  Chapter  3)  are  coupled,  which  is  characteristic  of  typical  FSI  problems. 


2.3  Response  to  Harmonic  Excitation 


Let  the  hannonic  excitation  with  circular  frequency  co  be  given  in  the  form 

{f{t)}=({fE}  +  {fp})exP{-i0)t)  (2.  i  3) 

where  |/£}  denotes  a  known  excitation  amplitude  which  may  be  acoustic  (say,  due  to  flow  noise) 

or  non-acoustic  in  origin,  and  j  fp  j  the  excitation  amplitude  associated  with  the  scattered  and 

radiated  acoustic  fields  for  elastic  scattering  problems.  Since  the  problem  under  study  is  linear,  the 
structural  response  to  harmonic  excitation  may  be  written  as 

=  {d}exp{-iat)  (2.14a) 

The  acoustic  pressure  field  in  the  fluid  domain  would  also  be  harmonic,  which  may  be  written, 
using  linear  superposition,  as 

P  (0  =  (Prs+  PRad )  exp  {-i  at)  (2. 14b) 


Here,  p^  denotes  the  total  complex  pressure  on  the  surface  of  the  body  resulting  from  rigid 

scattering  of  the  incident  field  by  the  body  located  in  the  waveguide,  which  of  course  is  obtained 
by  assuming  the  body  to  be  rigid.  This  involves  solving  a  rigid  scattering  problem  in  the 
waveguide,  which  is  discussed  in  Chapter  3.  Similarly,  pRad  denotes  the  complex  pressure  on  the 


body  surface  due  to  the  radiated  field,  which  is  caused  by  structural  motion  in  the  radial  direction 
only,  since  the  fluid  has  been  assumed  to  be  inviscid.  For  the  same  reason,  acoustic/hydrodynamic 
pressure  p  is  the  only  load  acting  on  the  structure  at  the  fluid-solid  interface  as  a  result  of  FSI.  So, 
in  view  of  Eq.  (2.9c),  the  virtual  work  done  by  these  pressure  forces  may  be  written  as 

(  »l+»2+»3 


swr 


I  Q«sd, 


exp  {-Hat) 


(2.15) 


V7=("i+«2+1)  ) 

where,  using  Eq.  (2.9b),  the  generalized  complex  pressure  force  may  be  written,  using  the  notation 
in  Eq.  (2.4c),  as 


j  =  («j  +n2  +l),(n,  +n2  +2), . (n]  +n2  +  nf) 

(2-16) 


=  0 


othewise 
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where  the  negative  sign  indicates  the  pressure  is  directed  opposite  to  the  outward  normal  of  the 
structural  surface.  Using  the  above,  the  harmonic  response  equation  for  the  structure  may  be 
written,  in  view  of  Eq.  (2.1 1)  as 


(-a>2[M}+[K}){d}  =  {Qr}  +  {QE}  (2.17) 

where  j  Qp  }  denotes  the  generalized  pressure  load  due  to  the  rigid  scattering  and  radiation  pressure 


acting  on  the  body  surface  (see  Eq.  (2.14b)),  which  will  be  dealt  with  in  Chapter  3,  and  {Qe}  the 
generalized  load  due  to  external  load  on  the  structure. 


2.4  Acoustic  Radiation  Loading 


Consider  the  radiated  pressure  pRad  in  Eq.  (2.16)  which  is  obtained  as  a  linear  combination  of 
complex  radiated  field  pRk  due  to  each  of  the  assumed  radial  modes,  say  Wt  in  Eq.  (2.4a).  This 
involves  solving  a  waveguide  propagation  problem,  in  which  the  body  is  located  in  the  waveguide 
and  the  normal  derivative  of  pressure  corresponding  to  A'-th  radial  mode  Wk  is  imposed  at  the  fluid- 
body  interface.  Note  that  this  condition  ensures  velocity  matching  at  the  interface.  The  above 
boundary  condition  for  the  acoustic  problem  may  be  written  in  general  as 
dpR,  d2wv 


=  ~P- 


(2.18) 


on  '  dt2 

where  n  denotes  the  normal  outward  of  the  fluid  domain,  p  the  density  of  the  fluid  and  the 
hannonic  radial  displacement  is  given  by 

wk  =  Wh  exp  (-i  cot)  (2. 19a) 

Then,  for  the  time-harmonic  case,  considering  the  k- th  radial  mode  mentioned  above,  in  view  of 
Eq.  (2.4a),  Eq.  (2.18)  reduces  as 

opRk 


dn 


■  =  pco-Wkdk 


(2.19b) 


dp 


It  is  convenient  to  solve  the  waveguide  problem  with  unit  amplitude  of  —  as  Wk  at  the  body- fluid 

dn 

interface.  Denoting  the  radiated  pressure  in  this  case  by  pRk0 ,  the  actual  radiated  field  due  to  the  k- 


th  mode  may  be  written,  including  the  modal  coefficient^ ,  as 

P  Rk  =  P®  P  Rkph  (2.20) 

Then,  in  view  of  Eq.  (2.16),  the  generalized  radiated  pressure  load,  with  all  the  assumed  modes 
included,  may  be  written  as 

«l+»2+»3  _ 

Qr,=\  X  pRkNjdS  j  =  nl+n2+l,nl  +  n2+2,...nl+n2+n3  (2.21a) 

$  k=n^+n2+l 


or 
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(2.21b) 


«i  +n2  +n  3 

e*=  I  Aa 

A=«j  +«2+l 

Added  mass  and  radiation  damping 

The  complex  loading  in  Eq.  (2.21)  can  be  decomposed  in  to  inertia  and  damping  load  components. 
The  inertia  load  component  may  be  written,  for  harmonic  case,  as 

-[j]  jjj  =  or  \A\{d}exp(-ia>t) 

where  the  negative  sign  is  based  on  D’Alembert’s  principle.  Here  [4]  has  the  unit  of  mass,  and 
hence  it  is  called  the  Added  Mass  matrix,  which  again  is  real  valued. 

The  damping  load  may  be  written,  for  the  hannonic  case,  as 
[i?]jjj  =  -ico[B]{d]  exp(-icot) 

where  [B]  has  the  unit  of  damping  constant,  and  hence  called  the  Radiation  Damping  Matrix, 
which  is  again  real  valued.  Note  that  both  [  A ]  and  [#]  are  symmetric  and  frequency  dependent. 


Thus,  the  complex  generalized  pressure  amplitude  in  Eq.  (2.21b)  may  be  written  as 

Fjk=co2[A{(o)\  -ieo\_B{co)\  (2.22) 

where, 

^=A'Re(^)  ;  (2.23a) 

CO ■  v  '  CO 

with  Re(-)  and  //«(■)  denoting  respectively  the  real  and  imaginary  parts.  In  actual  computation, 
the  factor  pco1  need  not  be  used  as  in  Eq.  (2.20),  and  hence  evaluate  the  integral  in  Eq.  (2.21a) 
using  pRk0  to  obtain^  .Then  Eq.  (2.23a)  may  be  written  as 

Ajk  =  PRe(Fjk )  ;  BJk  =  -pco lm ( F'jk )  (2.23b) 

The  primary  goal  of  waveguide  analysis,  with  the  radiating  body  included  in  the  fluid  domain 
(which  is  discussed  in  Chapter  3),  is  to  evaluate  added  mass  and  damping  parameters. 

In  view  of  Eqs.  (2.22)  and  (2.23),  the  harmonic  response  equation  in  (2.17)  may  be  written  as, 

(~co2  ([M]  +  p]) - ito[B]  +  [ *]) { d)  =  { Q '}  (2.24) 


where  [Q'}  denotes  generalized  loads  due  to  scattered  pressure  and  due  to  any  source  of  excitation, 

as  mentioned  under  Eq.  (2.13).  The  inviscid  fluid  model  used  does  not  obviously  give  rise  to  any 
tangential  fluid  dynamic  loading  on  the  body  fluid  interface  although  the  structural  motion 
considered  involves  both  radial  and  tangential  components  (see  section  2.2).  So  the  added  mass 
and  radiation  damping  matrices  in  Eq.  (2.24)  may  be  written  as 
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(2.25) 


0,  0 
0,  [A] 


0,  0 
0.  [B] 


where  the  trivial  entries  correspond  to  the  tangential  displacement  variables.  It  may  be  observed 
that  since  the  acoustic  radiation  loading  in  Eq.  (2.22)  depends  on  structural  motion,  its  effect  has 
been  appended  to  the  structural  mass  and  stiffness  terms  as  seen  in  Eq.  (2.24).  Solution  of  Eq. 
(2.24)  yields  the  structural  hannonic  motion  amplitude  array  [<f  j  associated  with  the  assumed 

modes  in  Eq.  (2.4a).  Having  obtained  the  structural  response  values,  the  radiated  field  in  the 
waveguide  can  readily  be  quantified  using  Eq.  (2.20). 


2.5  Conclusion 


The  structural  dynamic  model  discussed  in  Sec.  2.2  is  a  classical  fonn  which  is  especially  suited 
for  axisymmetric  structures,  and  it  results  in  a  compact  dynamic  model.  The  azimuthal/ 
circumferential  modes  would  obviously  be  Fourier  hannonics.  Legendre  functions  should  provide 
a  convenient  choice  for  the  meridional  modes.  It  might  be  noted  that  the  assumed  modes  approach 
discussed  above  would  readily  handle  mass  and  stiffness  variations  in  axisymmetric  shells,  as 
indicated  under  Eq.  (2.12).  If  it  is  desired  to  employ  FE  models  for  complex  axisymmetric  shell 
structures,  the  dry  modeshapes  may  first  be  derived  from  a  FE  model.  Then,  it  would  be  convenient 
to  obtain  a  least  squares  fit  to  these  meridional  modeshapes  using  assumed  modes  in  the  form  of 
Legendre  functions  and  then  employ  the  formulation  outlined  in  the  foregoing  sections. 
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Chapter  3 

Formulation  of  Acoustic  Scattering  and  Radiation  Problems 

in  3-D  Waveguides 


3.1  Introduction 

The  propagation  of  sound  waves  in  a  medium  may  be  analyzed  using  the  acoustic  fluid  model,  in 
which  the  fluid  is  considered  inviscid  but  compressible.  The  governing  equation  for  an  ocean 
acoustic  waveguide  (OAWG)  turns  out  to  be  the  wave  equation  in  terms  of  scalar  pressure, 
velocity  or  displacement  potential  (Jensen  et  ah,  2011).  The  wave  theoretical  model  can  handle 
heterogeneous  acoustic  properties,  backscatter  caused  by  obstructions  etc.  It  can  be  generalized  to 
include  attenuation  by  defining  a  complex  sound  speed.  Two  dimensional  propagation  in  a 
cylindrically  symmetric  waveguide  is  the  commonly  used  model  in  ocean  acoustics.  A  simplified 
model  for  three  dimensional  propagation  is  often  constructed  by  considering  such  2D  models  at  a 
set  of  discrete  azimuthal  angles  and  then  combining  the  various  2D  results.  Such  a  model  obviously 
does  not  account  for  azimuthal/horizontal  refraction/diffraction.  Three  dimensional  propagation 
arises  in  OAWG  due  to  factors  such  as  azimuthal  bathymetric  variations,  and  scattering  from 
objects  in  the  waveguide. 

In  this  chapter,  the  three-dimensional  (3-D)  Helmholtz  equation  governing  the  acoustic  field  in  a 
cylindrically  symmetric  waveguide  is  presented.  The  acoustic  field  is  expressed  in  a  Fourier  series 
in  the  circumferential/azimuthal  coordinates  and  the  governing  equation  reduced  to  a  2-D 
Helmholtz  equation  for  each  of  the  circumferential  harmonics.  The  axisymmetric  fluid  domain  is 
truncated  at  a  radial  boundary  in  the  far  field,  for  the  purpose  of  FE  modeling.  Hence,  special 
attention  is  focused  on  suitable  absorbing  boundary  conditions  on  the  truncation  boundary  instead 
of  using  the  classical  Sommerfield  radiation  condition,  which  is  a  limiting  condition. 

3.2  Mathematical  Model 

The  fluid  domain  il  defined  in  the  cylindrical  coordinate  system  (r,0,z\  constitutes  an  ocean 

acoustic  waveguide,  shown  schematically  in  Fig.  3.1,  where  z  denotes  the  vertical  coordinate,  r 
the  radial/range  co-ordinate  and  0  the  circumferential/azimuthal  coordinate. 
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The  waveguide,  unbounded  in  range  r,  is  bounded  at  the  top  by  the  plane  z  =  0 ,  which  is  the  air- 
sea  interface  or  free  surface  of  the  ocean (SF)  and  at  the  bottom  by  the  seabed  of  arbitrary 

topography  ( SB ) .  An  elastic  scatterer,  whose  surface  SIF  is  in  contact  with  the  fluid,  is  submerged 

in  the  waveguide.  The  wave  equation  for  the  heterogeneous  waveguide  may  be  written  in 
cylindrical  co-ordinates  as  (Jensen  et  ah,  2011), 


p(r,z,0)V. 


p(r,z,0) 


Vp(r,z,0,t ) 


+  - 


d2p(r,z,0,t ) 


(r,z,6) 


df 


=  0 


(3-1) 


where  p  denotes  the  time  dependent  acoustic  pressure,  p  the  density  of  the  medium  and  c  the 
sound  speed.  The  gradient  operator  in  cylindrical  coordinates  is  given  as 


/d_  }__d_  d\ 

\dr  ’  r  80  ’  dz] 


(3.2) 


Assuming  time-harmonic  motion  in  the  form p[r,z,0,t)=  p(r,z,0\e  l0" ,  where co  denotes  the 
circular  frequency  and  p  the  3-D  complex  acoustic  pressure,  the  3-D  Helmholtz  equation  follows 
from  Eq.  (3.1)  as 


p(r,z,0)V. 


p{r,z,0) 


Vp(r,z,0 ) 


+  - 


co 


( r,z,0 ) 


p(r,z,0)  =  0 


(3.3) 


Eq.  (3.3)  may  be  expanded  as 
\  Q  f  i 


P 


r  dr 


1  dp 
—  r  — 
\P  Sr 


+  P 


]_d_ 
r2  dO 


1  dp 
p  dO 


3  d 

+  P 


dz 


1  dp 
p  dz 


OJ  n 
+  —  p  =  0 

J  c_ 


(3.4) 


Let  the  solution  in  the  circumferential  coordinate  be  written  in  a  Fourier  series.  Thus,  the  total  field 
may  be  written  as, 


p{r,z,0)=  J  pm  (. r,z )  exp  (im6} 


(3.5) 
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Further,  it  is  assumed  that  p  and  c  are  independent  of  6.  This  implies  that  9  -  dependence  of  the 
acoustic  field  is  entirely  due  to  directional  sources,  scatterers,  and  radiators  which  are  located  in 
the  waveguide. 


It  is  important  to  recall  that  the  scatterer/radiator  has  been  assumed  to  be  axisymmetric.  Then,  in 
view  of  Eq.  (3.5),  Eq.  3.4  may  be  written  for  the  m -th  harmonic  in  a  decoupled  form  as 


1  r)  f 

p{r,z)~— 
r  or 


}_rdPn  A 
P  dr  J 


(  \  d 

+  P(r,z)— 

oz 


1  tyn 

\P  dz  J 


\  f 

+1 


CO " 


nr 


v 


\r,z) 


=  0 


(3.6) 


The  BVP  in  Eq.  (3.6)  should  satisfy  appropriate  B.C.  on  the  waveguide  boundaries  at  z  =  0  and 
z  =  D  .  A  prescribed  Neumann  B.C.  on  the  boundary  SN  may  be  written  as 

dPm  _ 


on 


on  S 


(3.7a) 


where  pmv  denotes  a  prescribed  nonnal  derivative  of  pressure.  The  rigid  seafloor  and  body-fluid 
interface  SIF ,  which  are  discussed  later  in  this  section,  would  fall  under  this  category. 


Alternatively,  a  prescribed  Dirichlet  B.C.  on  the  boundary  SD  may  be  written  as 

P,„  =  Pm<>  on  SD  (3.7b) 

For  pressure  release  surface  such  as  the  air-sea  interface  SF ,  pm0  =  0  . 

The  condition  on  an  impedance  boundary  Sr  may  be  written  as 

^  +  ap,=0onS,  (3.7c) 

on 

where  a  denotes  the  surface  impedance  coefficient. 

Apart  from  the  above  conditions,  the  solution  to  the  waveguide  problem  should  also  satisfy  a 
suitable  radiation  boundary  condition,  which  is  discussed  in  Sec.  4.3.  Eqs.  (3.6)  and  (3.7)  may  be 
adopted  to  study  a  variety  of  waveguide  scattering  problems,  as  discussed  below.  Consider  an 
axisymmetric  body  (with  a  vertical  axis  as  its  axis  of  revolution)  located  in  a  cylindrical  waveguide 
whose  reference  point  is  at  ( rB,zB  )  (see  Fig.  3. 1). 

The  acoustic  field  in  the  waveguide  due  to  a  point  source  will  be  scattered  by  this  body,  henceforth 
called  the  scatterer.  If  the  scatterer  is  deformable,  invoking  linear  superposition,  it  is  convenient  to 
consider  the  following  decomposition  (see  Eq.  (2.14b))  to  solve  this  problem: 

P  =  {Prs+ P  Rad)  (3-8a) 

where  pRs  denotes  the  complex  rigid  scattering  pressure,  pRad  the  complex  radiated  pressure 
owing  to  the  defonnability  of  body.  It  is  convenient  to  decompose  the  rigid  scattering  pressure  as 
a  sum  of  the  incident  pressure  field  pt  and  pressure  field  pRs0  due  scattered  waves  as, 

Pbs=(Pi+Psso)  (3-8b) 
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Note  that  the  decompositions  in  Eq.  (3.8)  lead  to  great  computational  convenience,  and  also 
provide  an  insight  into  elastic  scattering  phenomenon.  The  incident  pressure  in  Eq.  (3.8b)  indeed 
denotes  the  acoustic  field  caused  by  a  source  in  an  ‘empty  waveguide’,  i.e.,  the  waveguide  with 
the  scatterer  absent.  In  view  of  Eq.  (3.8),  the  following  waveguide  problems  can  be  identified: 

Case  1:  Empty  waveguide  (with  a  point  source) 


Considering  a  harmonic  omni-directional  point  source  at  r  =  0  and  z  =  zs  in  the  waveguide,  the 
forced  Helmholtz  equation  may  be  written,  using  Eq.  (3.6),  as  (Jensen  et  ah,  2011), 


p(r,z) 


2  q  f  i  A  a  (  i  a„  A  (  „2  A 


r  dr 


1  dp, 

— r  — — 
p  dr 


+  p(r,z) 


dz 


\_dp_i 

\P  dz  j 


+ 


CO 


[(r’z). 


_  _  c  Hr)S(Z~Zs) 

Pi  ~  ‘N 


2  nr 


(3.9) 


where  8  denotes  the  Dirac  delta  function.  Note  that  Eq.  (3.9)  governs  cylindrically  symmetric  wave 
propagation  in  a  range  and  depth  dependent  waveguide  containing  a  point  source  on  the  axis  of 
the  cylinder. 


Case  2:  Rigid  scattering  problem 


It  is  important  to  note  that  while  the  solution  to  Eq.  (3.9)  represents  an  axisymmetric  pressure  field, 
the  pressure  distribution  on  the  surface^  of  an  imaginary  scatterer  located  at  (rB,Zg) will  be 

nonaxisymmetric  with  respect  to  the  vertical  axis  of  revolution  of  the  scatterer.  Thus,  the  rigid 
scattering  problem,  despite  the  axisymmetry  of  the  scatterer,  is  indeed  a  nonaxisymmetric  problem. 
Fortunately,  choosing  a  suitable  coordinate  system,  the  incident  pressure  on  the  surface  S1F  may 
be  expressed  as  a  Fourier  series  as  in  Eq.  (3.5).  For  this  purpose,  consider  a  cylindrically  symmetric 
waveguide  in  a  new  coordinate  system,  (r,z)  with  the  axisymmetric  scatterer  located  at  7  =  0  and 
at  depth  zB( see  Fig.  3.1).  It  is  evident  that  the  radial  coordinate  in  Case  1  above  is  related  to  the 
radial  coordinate  here  as  r  =r-rB  .Then,  the  incident  pressure  field  may  be  written  following  Eq. 
(3.5)  as, 

m= oo 

Pj(r,z,6)=  Y  Pi.m(>7,z)exp(im0)  (3.10) 

m=— co 

And  the  pressure  due  scattered  waves  pRs0  in  Eq.  (3.8b)  may  be  written  as, 


PrsO  {r,Z,0)  =  Y  P  RsO,m  (r  >  z)  eXP  (  im0) 


(3.11) 


Then,  the  boundary  condition  on  the  surface  of  the  scatterer  is  given  by 

^PRsO,m  ^P\,m  ^  /-l 

— ^  =  — —on  SIF  (3.12) 

dn  dn 

where  n  is  normal  to  the  surface  of  the  scatterer.  The  governing  equation  for  this  waveguide 
problem  follows  from  Eq.  (3.6)  as, 
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/  \1  3 

p[r,z)~— 
r  or 


1  dpRs0n  A 
—  r - 

P  ^  j 


+  p{r,z)  — 

dz 


1  dp 


RsO.m 


p  dz 


+ 


CO 


m 


(r,z) 


P  RsO.m  ^ 


(3.13) 


The  forcing  for  this  problem  comes  from  the  interface  condition  in  Eq.  (3.12). 


Case  3:  Radiation  problem 


As  discussed  in  Sec.  2.3,  when  the  body  is  elastic,  it  vibrates  due  to  the  acoustic  pressure  loading 
on  it;  this  effect  can  be  represented  as  a  radiation  problem,  in  which  acoustic  waves  are  generated 
because  of  structural  motion  at  the  body-fluid  interface.  This  waveguide  also  has  7  as  the  range 
coordinate,  as  in  Case  2  above.  The  radiated  pressure  may  also  be  written  (see  Eq.  (3.11))  as, 

m= oo 

Pxad(7>z>0)  =  E  PRad,m(r,z)zxv{imO)  (3.14) 


Then,  the  boundary  condition  on  the  surface  of  the  radiator  is  given  by 

fy  rad , 


dn 


=  Pn 


(3.15) 


At  this  stage  it  is  important  to  recall  that  the  elastic  scattering  problem  being  a  coupled  one 
involving  fluid  dynamic  and  the  structural  dynamic  equations,  the  amplitude  of  normal  velocity 
field  can  assumed  to  be  unity,  but  its  distribution  on  the  surface  of  the  radiator  should  be  specified 
using  certain  assumed  modes.  Recalling  Sec.  2.5,  this  problem  can  be  tackled  following  one  of  the 
subsequent  approaches: 


1.  If  the  in-vacuo  structural  dynamic  problem  has  an  exact  solution,  the  mode  shapes  are 
known  analytically.  Several  problems  involving  the  vibration  of  cylindrical  and  spherical 
shells  are  discussed  in  Junger  and  Feit  (1972). 

2.  An  assumed  mode  approximation  may  be  used  when  the  structural  dynamic  equations  are 
not  amenable  to  exact  solution.  This  approach  has  been  outlined  in  Chapter  2. 


3.  For  complex  structural  forms,  the  dry  mode  shapes  of  vibration  may  be  generated  using  a 
finite  element  model.  Then,  in  order  to  follow  the  simplified  procedure  discussed  in 
Chapter  2  for  solving  the  radiation  problem,  an  assumed  approach  may  be  adopted.  The 
assumed  modes  may  be  fitted  to  the  mode  shapes  of  vibration  derived  from  the  FE  model 
employing  a  least  squares  method. 

The  governing  equation  for  this  waveguide  problem  follows  from  Eq.  (3.6)  as, 


/  n  1  d 
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\  f 

+  1 
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P  Rad.i 


=  0 


(3.16) 


The  forcing  for  this  problem  evidently  comes  from  the  interface  condition  in  Eq.  (3.15). 
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3.3  Radiation  Condition  at  the  Truncation  Boundary 

Wave  propagation  problems  require  suitable  radiation  boundary  condition  to  ensure  that  the  waves 
are  truly  outgoing.  In  linear  problems,  the  well-known  Sommer  field  radiation  condition,  which  is 
a  limiting  condition,  is  commonly  used.  Domain  discretization  methods  such  as  the  finite 
difference  and  FE  methods  applied  to  wave  propagation  problems  necessarily  require  appropriate 
truncation  of  the  problem  domain.  In  order  to  approximate  the  radiation  B.C  on  such  truncation 
boundaries,  several  absorbing  or  nonreflecting  boundary  condition  have  been  developed.  In  the 
present  study,  following  Fix  and  Marin  (1978)  and  Vendhan  et.  al  (2010),  the  boundary  damper 
operators  due  to  Bayliss,  Gunzberger  and  Turkel  (BGT)  (1982)  are  adopted  for  use  in  the  FE 
model. 


The  BGT  damper  equation  for  a  cylindrically  symmetric  wave  may  be  written  in  a  fonn  similar  to 
Eq.  (3.7c).  The  first  order  cylindrical  damper  equation  may  be  written  on  the  truncation  boundary 
SR  as 


^L  +  apm  =0onSR  (3.17) 

or 

where  for  exterior  hannonic  acoustic  problems,  a  =  -ik ,  the  medium  wavenumber.  On  the  other 
hand,  in  waveguide  acoustics,  a  finite  number  (say  n)  of  radial  modes,  with  associated  wave 
numbers,  travel  outward,  and  hence  Eq.  (3.17)  needs  to  be  satisfied  for  each  of  them  (Fix  and 
Marin,  1978  and  Vendhan  et  al.,  2010).  Thus,  for  cylindrical  waveguides,  Eq.  (3.17)  may  be 
written  as 


dPm] 

dr 


+  a,Pn*  =  0  011  SR 


(3.18a) 


where 


aj  = 


r--ik 
2  r  1 


>7  =  1,2,3...  (3.18b) 

\^>  j 

In  Eq.  (3.18b),  k.,j  =  1,2,3... ,  n  denote  the  radial/horizontal  wavenumbers.  For  depth  dependent 


waveguides,  the  horizontal  wavenumbers  krj  and  the  associated  depth  modes  /  (z)  can  in  general 

be  obtained  numerically.  It  is  well  known  the  nonnal  mode  solution  for  such  a  waveguide  may  be 
written,  at  a  given  range,  as 

Pm(z)  =  iajfj(z )  (3-19) 

j= i 

where  the  modal  coefficients  a j  may  be  obtained  using  the  orthogonality  of  the  depth 

modes/eigenmodes.  In  view  of  Eq.  (3.19),  the  absorbing  boundary  conditions  in  Eq.  (3.18a)  may 
be  written  as 


dP,n 

dr 


+  X  */*/;(*)  =  ° 


7=1 


(3.20) 
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3.4  Variational  Formulation 


In  order  to  facilitate  finite  element  modeling,  a  variational  fonnulation  equivalent  to  Eq.  (3.6) 
together  with  appropriate  boundary  conditions  is  chosen  here.  Consider  the  functional  given  by 


l(r,z )  =  —  [ — - — 
V  1  2 1  P(r,z) 


+r 


dr 


dp ^ 
dz  j 


rdr  dz  - 


iJ— 
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\r,z) 


p2  rdr  dz 
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SD+Sp, 


_ dp__  _ 

p(r,z)dnP  2  j=l  p(r,z) 


|  (a  .p2^  rdz 


(3.21) 


For  Case  1  defined  Sec.  3.2, 

m  =  0  and  p  =  pIm  (3.22a) 

where  the  point  source  at  r  =  0  provides  the  excitation. 

For  Case  2  defined  in  Sec.  3.2, 

P  =  Prs  o,m  (3.22b) 

the  /;z-th  harmonic  component  of  scattered  pressure.  The  B.C.  on  the  scattering  body  is  given  by 
Eq.  (3.12). 

For  Case  3  defined  in  Sec.  3.2, 


P  P Rad, m 


(3.22c) 


the  /H-th  radiated  acoustic  pressure.  The  body  B.C.  given  by  Eq.  (3.15). 

It  can  be  shown  that  stationarity  condition 

<57  =  0  (3.23) 

For  the  functional  in  Eq.  (3.21)  leads  to  the  governing  equation  and  the  associated  boundary 
conditions  discussed  above. 


3.5  Steps  in  Elastic  Scattering  Analysis 

In  order  to  provide  a  better  understanding  of  the  computational  model  discussed  above,  a  flow 
chart  of  the  various  steps  is  given  in  Figs.  3.2a  and  3.2b. 
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STEP  1 

Waveguide  coordinates:  (r,z) 

Solve  empty  waveguide  with  point  source  at  r  =  0 
to  find  Pj 

(This  is  an  axisymmetric  WG  analysis) 

dp 

Find  — -  on  the  imaginary  surface  of  the  scatterer 
dn 

(This  will  be  nonaxisymmetric  w.r.t.  body  axis)) 


STEP  2 

Waveguide  coordinates:  (r,z) 

Body  centered  at  r  =  0 

Perform  nonaxisymmetric  WG  analysis  to  find 
rigid  scattered  pressure  pRS  . 

Body  boundary  condition  (Eq.  (3.12)) 

dpRS(r,z)_  dpI  (r,z) 
dn  dn 


Find  generalized  load  vector  due  to  scattered 
pressure  (see  Eq.  2.16) 

Fig.  3.2a  Steps  in  Rigid  Scattering  Analysis 
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STEP  3  (repeated  for  m  =  1,2,..  ) 

Waveguide  coordinates  (r,z) 

Body  centered  at  r  =  0 

Perform  nonaxisymmetric  WG  analysis  to  find 
radiated  pressure  pRadJ  due  to  assumed  mode  j  with 

unit  amplitude  of  pressure  gradient  on  the  scatterer 
surface. 

Body  boundary  condition  :  ^ ^  ^  -  f  (r,z) 

dn 


Find  acoustic  radiation  load  matrix  FJk  (see  Eq. 

2.21) 

Derive  added  mass  and  radiation  damping  matrices 
from  Fjk  (see  Eq.  2.22) 


Fig.  3.2b  Steps  for  Setting  up  Added  Mass  and 
Radiation  Damping  Matrices _ 
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Chapter  4 

Finite  Element  Analysis  of  Depth  and  Range  Dependent  Waveguides 

4.1  Finite  Element  Model 

The  variational  problem  in  Eq.  (3.21)  may  be  solved  by  extending  the  finite  element  model  in 
Vendhan  et  al  (2010)  to  nonzero  circumferential  harmonics,  i.e.,  m  >  0  .  An  8-node  isoperimetric 
quadrilateral  element  has  been  adopted  in  the  FE  fonnulation.  The  acoustic  pressure  p  in  a  finite 
element  may  be  written  as 

m 

p(r,  z)  =  X  PejNj  (£  £)  =  [#]{&}  (4.1) 

7=1 

where  N.  denotes  a  polynomial  shape  function  in  terms  of  the  parametric  coordinates  (£,$)  in  the 
(r,z)plane  and  pej  the  acoustic  pressure  at  the  j- th  node  of  an  m-node  finite  element.  For  details  of 

isoparametric  quadrilateral  element  formulation  see  Cook  et  al.  (2004).  Upon  using  Eq.  (4.1)  in 
Eq.  (3.21),  the  following  discrete  approximation  for  the  functional  is  obtained: 

I(Pe)a\{PeY  i[Ke]-[Me\){Pe}  +  \Y.{Pem}T\^em'\{Pern}-{Pe}T  {fe}  (4'2) 

The  stationarity  condition  (see  Eq.  (3.23))  on  the  potential  I (pe) in  Eq.  (4.2)  should  be  subjected 
to  the  constraint  so  as  to  impose  the  boundary  damper  condition  in  Eq.  (3.18).  The  damper  equation 
in  (3 . 1 8)  may  be  written  (at  r  =  rB  ,  the  radiation  boundary  SR  )  as 


[C]  ( p(rB,z ),  a,  a2,  ...an)T  =  0  (4.3a) 

where, 

[c]  =  [1,  -/  (Z),  -/j  (z).  ■•■-/.  (z)]  (4.3b) 

In  Eq.  4.3a,  <5^  are  the  unknown  modal  coefficients  associated  with  the  depth  modes /y(z) . 


Adopting  the  penalty  function  approach  discussed  in  Vendhan  et  al.  (2010),  this  may  be  achieved 
by  constructing  a  modified  potential  as 

’'  =  ‘  +  {{p'.}T[C^[Pr][C,]{p'.}  <4-4a> 

i.e. 


where  the  expanded  nodal  dof  array  is  given  by 

{*H{a 


(4.4b) 


(4.5) 
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The  various  matrices  in  Eq.  (4.4b)  are  given  by 


’[*.]  o' 

O  0 


■\K}  = 


\M,\  o' 
0  0 


0  0 

0  ft] 


;  (/,'}  = 


0 


[JC,]=J-V[JV]rV[AT]dn 

n.  P 


[Me]=\-[Nf[N]dn 
Q(.  P 

[Re]  =  diag[Rel,Re2,...ReM  ] 

K.  =  {/_}r[s..]{/-}.R»]  =  a»  i  -[Vf[Jv]<is 


P 


{fzm  }  =  (fm  (Z1 )  >  -L  (Z2  )  > . /„  (Z„  ))' 


P 


(4.6a) 

(4.6b) 

(4.6c) 

(4.6d) 

(4.6e) 

(4.6f) 

(4.6g) 


where  [Ke],  [Me]  and  [ Re ]  are  fonnally  known  as  the  element  stiffness,  mass  and  radiation 
damping  matrices,  and  [fe]  denotes  the  element  load  vector.  The  stationarity  condition  in  Eq. 
(3.23)  is  replaced  by  the  following  condition 


dl' 

s{pI 


=  0 


(4.7) 


which  leads  to  the  matrix  equation 

(K]-[M;]+[^]+[c,f[/?,][c,]){K}={/:}  (4.8) 

The  finite  element  assemblage  operation  using  the  element  equations  in  (4.6)  leads  to  the  global 
FE  equation  given  by 

(ft]  -ft ']  +  ft]  +  [C']){/>'}  =  {/'}  (4.9) 

where 

ftftlft]  ;ft']=Zft:hft]=Ift]  (4.ioa) 

e  e  e 

[C']  =  Z[C.]r[/?„][C.] ;  {/'}  =  S{/,'}  (4.10b) 

e  e 

with  ^  denoting  the  assemblage  operation. 
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4.2  Solution  of  the  FE  Equations 


The  linear  matrix  equation  in  (4.9)  resulting  from  the  FE  model  for  waveguide  propagation  may 
be  written  as 

MW={/}  (4-11) 

The  matrix  [A]  is  symmetric  and  positive  definite.  A  FORTRAN  90  code  has  been  developed  as  a 
part  of  the  present  study,  implementing  the  FE  model  for  the  modal  Helmholtz  equation  for  solving 
both  rigid  scattering  and  radiation  problems.  A  block-storage  banded  Gauss  solver  has  been 
employed  to  solve  Eq.  (4.11).  For  large  FE  models,  iterative  solvers  based  on  the  conjugate 
gradient  method  would  be  computationally  more  efficient.  However,  it  should  be  remembered  that 
the  capability  of  the  Gauss  solver  to  handle  multiple  load  vectors  with  a  single  factorization  could 
be  exploited  for  solving  radiation  problems  involving  many  assumed  modes.  As  indicated  in  Eq. 
(3.12),  the  normal  derivative  of  the  incident  pressure  field  needs  to  be  calculated  from  the  FE 
solution.  Since  the  FE  model  in  Eq.  (4.1)  denotes  a  C0  model,  calculation  of  derivatives  requires 

special  procedures  which  are  analogous  to  methods  employed  for  stress  calculation  from  FE 
solution  to  elasticity  problems  (Cook  et  ah,  2004). 


4.3  Depth  Eigenproblem 


The  mathematical  model  that  fonns  the  basis  for  the  FE  model  is  based  on  the  assumption  that  in 
the  vicinity  of  the  truncation  boundary  and  beyond  the  waveguide  has  constant  depth  and  the 
acoustic  properties  are  only  depth  dependent.  Thus  the  FE  model  requires  depth  eigensolution  at 
the  truncation  boundary.  For  this  purpose,  Eq.  (3.6)  may  first  be  reduced  to  the  following  form 
assuming  depth  dependent  properties: 


?(>■)-- 

'  r  dr 


I  r5Pn 


+  p{r)  — 
v  ’  dz 


1  tym 

yP  dz 


^  ^  (Q2 


P  Sr 

Assuming  a  variable  separable  solution  in  the  form 

Pm{r,z)  =  R(r)z(z) 

Eq.  (4.12)  may  be  written  as 


+ 


nr 


c2(z)  r2 


p  =0 

i  m 


P 


ld_ 

r  dr 


1  &R  'i 

—  r  — 

P  Sr  J 


nr  d 

— R  P 


dz 


1  dZ  )  a2 


p  dz 


R 


=  0 


(4.12) 


(4.13) 


(4.14) 
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Then,  the  depth  eigenproblem  for  the  axisymmetric  waveguide  with  nonaxisymmetric  acoustic 
field  in  the  form  given  by  Eq.  (4.13),  of  course  with  only  depth  dependence,  may  be  deduced  from 
Eq.  (4.14)  as 


P 


dz 


f  1  dZ ^ 

— 

+ 

—r-K 

Ip  & ) 

lA2  rJ 

Z  =  0 


(4.15) 


where  ^.denotes  the  radial  wavenumber.  Eq.  (4.15)  is  evidently  identical  to  the  case  of 
axisymmetric  wave  field.  The  depth  modes  Z.(  =  )  and  radial  wavenumbers  krJ  may  readily  be 

obtained  using  the  KRAKEN  code  [Porter,  2001],  which  employs  a  finite  difference  model  to 
discretize  the  depth  equation  in  (4.15).  However,  in  the  present  study,  a  Rayleigh-Ritz  (RR)  model 
has  been  developed  to  obtain  the  eigensolution,  in  which  the  depth  modes  appear  in  a  compact 
form  as  a  linear  combination  of  the  isovelocity  mode  shapes.  Such  a  form  for  the  mode  shapes  has 
been  found  to  be  more  convenient  to  use  in  the  proposed  FE  model.  The  details  of  the  RR  model 
for  the  depth  eigenvalue  problem  and  numerical  validation  are  presented  in  Chapter  5. 
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Chapter  5 


Rayleigh-Ritz  Approximation  to  Depth  Eigenvalue  Problem  of 
Shallow  Water  Acoustic  Waveguides 


5.1  Introduction 

The  sound  speed  in  an  ocean  acoustic  waveguide  is  in  general  both  depth  and  range  dependent. 
Depth  dependence  is  considered  very  important  because  it  is  responsible  for  many  interesting 
phenomena  in  waveguide  propagation.  The  two  well-known  methods  that  have  been  developed 
to  study  acoustic  waves  in  depth  dependent  waveguides  are  the  fast- field  technique  and  the  normal 
mode  expansion  (Frisk,  1994;  Jensen  et  ah,  2011),  the  latter  being  the  subject  matter  of  the  present 
paper.  The  normal  mode  approach  consists  of  first  solving  the  depth  eigenvalue  problem  for  a 
given  sound  speed  profile  to  obtain  the  radial  wavenumbers  and  the  associated  depth  modes,  which 
respectively  are  the  eigenvalues  and  eigenfunctions.  The  depth  eigenproblem  could  be  solved 
exactly  only  for  a  few  cases.  So,  several  numerical  methods  have  been  developed  to  solve  the 
depth  problem  (Jensen  et  ah,  2011).  The  KRAKEN  code  (Porter,  2001)  employs  a  finite  difference 
model  for  the  depth  equation  and  the  resulting  algebraic  eigenproblem  has  been  solved  using  a 
combination  of  iterative  techniques  and  Richardson  extrapolation  to  obtain  the  radial 
wavenumbers  and  modal  vectors  to  a  great  degree  of  precision  (Porter  and  Reiss,  1984).  Given  the 
depth  eigensolution,  the  acoustic  field  in  the  waveguide  can  be  obtained  using  the  well-known 
nonnal  mode  approach  (Kinsler  et  ah,  2000,  Porter,  2001;  Jensen,  2011).  The  KRAKEN  code  has 
been  used  extensively  in  the  ocean  acoustic  community  for  both  shallow  and  deep  water 
waveguides.  In  the  finite  element  model  for  depth  and  range  dependent  waveguides  (Fix  and 
Marin,  1978;  Vendhan  et  al.,  2010),  the  eigenvalue  solution  of  the  depth  problem  is  required  for 
imposition  of  the  radiation  condition  at  the  truncation  boundary.  In  the  Couple  method  (Evans, 
1984),  the  depth  eigenvalue  solution  is  required  at  every  range  independent  strip,  a  series  of  which 
makeup  a  range  dependent  waveguide.  For  the  foregoing  applications,  it  would  be  convenient  to 
have  the  depth  modes  in  a  compact  form.  The  main  thrust  of  the  present  paper  is  to  explore  this 
aspect  with  specific  reference  to  shallow  water  waveguides. 

The  depth  eigenvalue  problem  can  be  cast  in  a  variational  fonn  by  suitably  defining  a  functional. 
Then,  the  classical  Rayleigh-Ritz  (RR)  method,  widely  used  to  find  approximate  solutions  of 
eigenvalue  problems  in  structural  mechanics,  may  be  employed  to  find  a  variational  approximation 
to  the  eigenvalue  solution  of  the  depth  problem  in  ocean  acoustic  waveguides.  The  depth  modes 
obtained  would  have  a  more  compact  form  than  those  derived  using  finite  difference  or  even  finite 
element  methods.  There  seems  to  be  no  study  dealing  with  this  approach.  The  present  paper 
provides  a  RR  model  for  the  depth  eigenvalue  problem  and  demonstrates  its  accuracy  for  shallow 
water  waveguides.  The  numerical  features  of  the  RR  model  are  brought  out  by  applying  it  to  both 
single  and  two  layer  waveguides. 
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The  depth  equation  for  an  acoustic  waveguide  overlying  a  fluid  halfspace  is  presented  in  Sec.  5.2. 
A  variational  fonnulation  for  the  depth  problem  is  presented  and  a  RR  approximation  to  the  depth 
eigensolution  is  discussed  in  Sec.  5.3.  Sec.  5.4  provides  a  brief  overview  of  the  numerical  aspects 
of  the  RR  model,  and  its  application  to  both  single  and  two  layer  shallow  water  waveguides.  Salient 
features  of  the  RR  model  are  discussed  and  results  based  on  this  model  are  compared  with  the 
exact  solution  and  numerical  results  from  KRAKEN.  Sec.  5.5  reviews  the  RR  model  presented 
and  its  possible  extensions. 

5.2  Mathematical  Model 


Consider  a  cylindrically  symmetric  waveguide  having  depth  dependent  density  p  and  sound  speed 
c  .  A  homogeneous  acoustic  fluid  half-space  is  assumed  to  underlie  the  OAWG.  The  OAWG  is 
assumed  to  be  the  inhomogeneous,  and  the  pseudo  Helmholtz  equation  governing  the  linear 
harmonic  acoustic  pressure  field  p(r,z)  in  the  waveguide.  It  is  given  in  cylindrical  coordinates 
(r,z)  as  (Frisk,  1994;  Jensen  et  al.,  2011) 


]_8_ 

r  dr 


dp 
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+  p{z) 
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dz 


1  dp 
p{z)  dz 


+  - 
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C  (z) 


8{r)8{z-zs) 
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(5.1) 


where  r  denotes  the  range  coordinate  and  z  the  depth  coordinate  as  shown  in  Fig.  5 . 1 ,  and  the  r.  h.s. 
denotes  a  point  source  of  unit  amplitude  located  at  r  =  0  and  z  =  zs,  with  8  denoting  the  Dirac 
delta  function. 
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Figure  5.1  A  two  layer  waveguide 

Eq.  (5.1)  can  also  be  applied  to  problems  with  attenuation  by  introducing  a  complex  sound  speed. 
In  order  to  set  up  the  nonnal  mode  solution  of  Eq.  (5.1),  a  variable  separable  solution  of  the 
homogeneous  form  of  Eq.  (5.1)  may  be  written  as 

p(r,z)  =  R(r)Z(z)  (5.2) 
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Then,  upon  using  Eq.  (5.2)  in  the  homogeneous  form  of  Eq.  (5.1),  the  following  ordinary 
differential  equations  are  obtained  (see  Sec.  4.3): 


1  d 
r  dr 

PO) 


dR(r) 

dr 


+  k2R(r)  =  0 


d_ 

dz 


1  dZ{z) 
p{z)  dz 


+ 


f  CD2 


-k: 


Z(z)  =  0 


(5.3) 


(5.4) 


Vc"(z) 

where  co  denotes  the  circular  frequency  and  k]  the  separation  constant,  which  turns  out  to  be  the 
square  of  the  radial/horizontal  wavenumber.  Eq.  (5.3)  evidently  pertains  to  the  radial/horizontal 
modesR(r)  and  Eq.  (5.4),  which  is  of  main  concern  in  the  present  paper,  pertains  to  the  depth 

modes  Z  (z) .  Choosing  a  pressure  release  boundary  at  the  top  ( z  —  0  )  and  a  mixed/Robin  boundary 
condition  at  the  seabed  (z  =  Dx ),  they  may  be  written  as  (Jensen,  2011;  Porter,  2001), 

Z(0)  =  0  (5.5) 


Z{Dx)  +  ^l—  =  0 

p  dz 


(5.6a) 


where 


g 


(<*;)  =  pJ. 


^ r  2 

V  c  7 


(5.6b) 


with  ph  denoting  the  density  of  the  acoustic  fluid  in  the  isovelocity  half-space  underlying  the  water 

column.  Eq.  (5.6)  facilitates  replacing  the  half-space  in  the  Pekeris  waveguide  (Pekeris,  1948)  by 
means  of  an  impedance  type  boundary  condition.  It  is  interesting  to  note  that  Eq.  (5.4)  together 
with  the  homogeneous  boundary  conditions  in  Eq.  (5.5)  and  Eq.  (5.6a)  does  not  constitute  a  proper 
Sturm-Liouville  problem  because  Eq.  (5.6a)  depends  on  the  unknown  eigenvalue  kr2 .  In  the 
KRAKEN  code,  a  finite  difference  model  is  employed  to  solve  Eq.  (5.4)  together  with  the 
boundary  conditions  in  Eq.  (5.5)  and  Eq.  (5.6). 


As  an  alternative  to  Eq.  (5.6a),  the  depth  equation  for  the  fluid  half-space  is  considered  here. 
Noting  Eq.  (5.1),  for  the  half-space, 


P,XZ) 


d_ 

dz 


1  dZb{z) 
ph{z)  dz 


+ 


co~ 


\ciXz) 


-k: 


Zb(z)  =  0,  Dx<z<  CO 


(5.7) 


where  Zh  denotes  the  depth  function  in  the  fluid  half-space  having  depth  dependent  density  ph  and 
sound  speed  ch  .The  interface  conditions  at  the  seabed  are  given  by  the  kinematic  and  dynamic 
conditions, 
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Z(Dl)  =  Zh(Dl) 


(5.8a) 


1  dZ(D{)  _  1  dZh(Dl)  ^  8bN 

p  dz  pb  dz 

In  addition,  the  depth  mode  Zh  should  remain  bounded  as  z  — »  <x> .  The  objective  here  is  to  consider 
a  variational  formulation  for  Eq.  (5.4)  and  Eq.  (5.7)  together  with  appropriate  boundary  conditions 
and  obtain  a  RR  approximation  to  the  depth  problem. 

5.3  Variational  Formulation  and  Rayleigh-Ritz  Approximation 

A  variational  formulation  equivalent  to  the  boundary  value  problem  in  Sec.  5.2  is  sought  now.  The 
operator  being  symmetric,  there  exists  a  functional,  the  variation  of  which  leads  to  Eq.  (5.4)  and 
appropriate  boundary  conditions,  and  similarly  for  the  half-space.  Consider  the  functional  Y\(Z) 

and  Ylb(Zh)  defined  respectively  in  the  water  column  and  the  half-space  as 


l  D' 

n(zKl 


p(z) 


dZ_ 
ydz  j 


(D~ 


1 


,  -Z2+ - k;(z)Z2 

p(z)c  (z)  p(z ) 


\dz- 


p(z) 


Z.Z 


(5.9) 


A 

n,=-| 


D2 

f 

1  | 

(dZt\ 

j 

A 

_a(zV 

y  dz  y 

co 


zf,+- 


2ry2 


Pb(z)c~(z )  ph(z) 


-ktZ 


dz 


a,(a>  2  pb 


(5.10) 


1  „Z2b(D2 ) 


,D2-+co 


At  the  interface  z  =  Dx  between  the  water  column  and  the  half-space,  the  conditions  noted  in  Sec. 
5.2  must  be  imposed.  In  view  of  Eq.  (5.8b),  this  can  be  achieved  by  setting  in  Eq.  (5.9) 

=  (5.11) 

P  Pb  dl 


and  in  Eq.  (5.10), 


Pb 


1  dZ(D{) 
p  dz 


(5.12) 


In  addition,  Eq.  (5.8a)  should  be  imposed.  Then,  it  can  easily  be  shown  that  the  variational 
condition  <7 FI  =  0  leads  to  Eq.  (5.4)  and  the  boundary  conditions  in  Eq.  (5.5)  as  well  as  the  interface 
condition  in  Eq.  (5.11),  where  £  denotes  the  first  variation.  Similarly,  the  variational  condition 
SUb=0  leads  to  Eq.  (5.7),  and  the  interface  conditions  in  Eq.  (5.8a)  and  Eq.  (5.12).  In  addition, 
at  z  =  D2  we  obtain  the  condition 
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^j^+f3Zb(D2)=  0 

dz 


(5.13a) 


where,  noting  Eq.  (5.6),  for  a  homogeneous  half-space, 


(5.13b) 


Note  that  for  waveguides  with  a  finite  depth  fluid  layer  under  the  water  column,  the  last  term  in 
Eq.  (5.10)  is  deleted  and  D2  would  be  finite. 


We  now  seek  an  assumed  mode  solution  with  n  tenns  to  the  above  variational  problem  in  the  form 
Z{z)»  j^^.(z),  0  <z<Dx  (5.14a) 

M 

Zb(z)  ~  (5.14b) 

7=1 

where  ^  denotes  a  known  coordinate  function  satisfying  the  kinematic  boundary  condition  in  the 
water  column  and  (/>,  an  unknown  constant,  and  their  counterparts  with  suffix  b  correspond  to  the 

half-space.  Let  the  two  sets  of  coordinate  function  above  be  such  that  they  satisfy  the  interface 
conditions  in  Eq.  (5.8)  and  the  condition  in  Eq.  (5.13).  Such  functions  may  readily  be  constructed 
by  solving  a  two  layer  depth  problem,  (which  is  nothing  but  the  Pekeris  waveguide),  with  an 
appropriate  choice  of  constant  velocity  and  density  in  the  water  column  and  the  half-space.  This 
approach  has  been  adopted  here.  Then,  it  follows  that  the  continuity  of  pressure  field  at  the 
interface  z  =  Dx  implies  that  in  Eq.  (5.14) 

(j>j  -(j)j ,  n  =  n1=n2  (5.15) 

Then,  using  Eq.  (5.14)  in  the  functional  in  Eq.  (5.9)  and  Eq.  (5.10)  and  combining  them,  an 
algebraic  approximation  for  the  functional  is  obtained  as 

n=n+nl=t{?}T[A-']{f}-io!{f}r[jf"]{f}+lk;{f}T[jf'"]{?}  (5.i6> 

where 
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1  dy/i  dy/ j 


yO(z)  Jz  Jz 

1 


2 

Jz  +  J 


1  dYbi  dVbj 


o  P(z)c  (z) 


y/^jdz  +  J 


pb  (z)  Jz  Jz 
■» 
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Jz 


Pb(z)cl(z) 


VbWhA2 


A  j  A  j 

a = i  + j 


A  oo 


(5.17) 


It  has  been  assumed  in  the  above  the  contribution  due  to  the  last  two  tenns  in  Eq.  5. 10  is  negligible 
as  D2  — »  oo  . 


The  variational  condition  is  now  replaced  by  the  condition 


on 

dcj)j 


=  0, 


j  =  1,2,..., ft 


Eq.  (5.18)  yields  a  symmetric  algebraic  eigenproblem  given  by 


nr 


(5.18) 


(5.19) 


The  eigensolution  of  Eq.  (5.19)  may  be  denoted  as 

(^,{^})j  =  l,2,..n  (5.20) 

Having  obtained  the  eigenvalues  k2rj  and  the  eigenvectors  ,  the  eigenfunctions/depth  modes 
may  be  written,  using  Eq.  (5.14),  as 

r(z)={?°’f  h(z)}  (5.2i) 


where 


{y/(z)}T  =<^(z),^2(z),..^„(z)>  (5.22) 

Eq.  (5.21)  provides  a  compact  for  the  depth  modes  that  are  convenient  to  employ  in  FE  models 
such  as  those  in  Fix  and  Marin  (1978),  Jensen  et  al.  (2011),  and  Vendhan  et  al.  (2010)  for 
approximating  the  radiation  condition  at  the  truncation  boundary.  The  depth  modes  obtained  can 
of  course  be  used  to  set  up  the  nonnal  mode  solution  to  the  forced  Helmholtz  equation  in  Eq.  (5.1). 
Since  the  eigenvectors  in  Eq.  (5.19)  are  -orthogonal,  it  can  easily  be  shown  the 

eigenfunctions  in  Eq.  (5.21)  satisfy  the  following  orthogonality  condition: 
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1 


(5.23) 


A 

I 


a(z) 


Z;Z;<iz  =  0 


In  terms  of  finite  element  tenninology,  the  RR  model  for  each  layer  may  be  looked  upon  as  a  super 
element  with  C,  continuity  at  the  inter-element  boundary  and  the  operation  leading  to  Eq.  (5.16) 
is  equivalent  to  element  assemblage  operation. 


5.4  Numerical  Analysis  and  Discussion 

The  Rayleigh-Ritz  model  in  Sec.  5.3  has  been  implemented  in  MATLAB.  The  first  task  is  to 
compute  the  symmetric  matrices  [ K1  ] ,  [ATW]  and  ^K1"  J  in  Eq.  (5.17).  The  associated  integrals 

have  been  evaluated  using  composite  Gauss  quadrature  functions  available  in  MATLAB.  The  next 
task  is  to  find  the  eigensolution  to  Eq.  (5.19).  For  problems  with  no  attenuation,  the  real 
eigenvalues  have  been  obtained  employing  the  bisection  method.  For  problems  with  attenuation, 
approximations  to  the  complex  roots  have  been  obtained  using  a  search  procedure  (Buckingham 
and  Giddens,  2006)  and  the  eigenvalues  refined  employing  Newton-Raphson  iteration.  In  all  cases, 
the  eigenvectors  are  obtained  using  inverse  iteration. 

The  Rayleigh-Ritz  model  in  Sec.  5.3  has  been  applied  first  to  single  layer  isovelocity  waveguide 
examples  without  attenuation  for  which  exact  solutions  are  available.  The  eigenvalues  computed 
(not  reproduced  here  for  brevity)  compared  with  the  exact  results  better  than  8  or  9  significant 
digits.  This  in  essence  provides  validation  of  the  numerical  accuracy  of  the  matrices  Kl  ] ,  \^Kn  J 

and  ^//y]in  Eq.  (5.17),  evaluated  using  numerical  quadrature  in  MATLAB,  especially  for  the 

half-space.  Then,  three  different  sound  speed  profiles  (Sazontov;  2002;  Schneider,  1977)  shown 
in  Fig.  5.2a-c  have  been  chosen  to  evaluate  the  accuracy  of  the  RR  model.  Attenuation  in  the  fluid 
half-space  has  also  been  considered.  Different  sets  of  RR  approximations  have  been  obtained  by 
varying  the  number  of  assumed  modes  n  in  Eq.  (5.15).  The  results  for  n  =  2 nP,  where  nP  denotes 
the  number  of  propagating  modes  turned  out  to  be  of  good  accuracy  and  hence  reproduced  in 
Tables  1-6.  All  the  examples  studied  have  also  been  solved  using  a  FORTRAN  version  of  the 
KRAKEN  code  (Acoustics  Tool  Box),  and  the  results  compared  with  the  present  RR  model.  It  is 
seen  that  most  results  compare  up  to  5  significant  digits,  and  often  better,  which  is  deemed  very 
good.  Even  the  imaginary  parts  of  the  wavenumbers  for  the  attenuating  half-space  problem 
compare  fairly  well.  The  modeshapes  obtained  by  the  two  methods  are  indistinguishable  to  plotting 
accuracy,  and  hence  not  reproduced  here. 

The  elements  of  the  RR  matrix,  both  in  the  finite  layer  and  the  semi-infinite  layer,  could  be 
computed  to  a  high  degree  of  precision  using  MATLAB  tools.  The  excellent  perfonnance  of  the 
RR  model  for  the  depth  eigenproblem  may  be  attributed  to  the  following: 
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a)  The  modeshapes  of  an  isovelocity  waveguide  have  been  chosen  as  trial  functions, 
which  satisfy  appropriate  interface  conditions  at  the  sea  floor,  and  the  boundary 
conditions  at  the  free  surface  and  finite  bottom.  This  renders  the  RR  matrix  highly 
diagonally  dominant,  which  also  greatly  aids  in  numerical  evaluation  of  the 
eigensolution. 

b)  For  ocean  waveguides,  the  depth  variation  of  the  sound  speed  is  normally  only  a  small 
percentage  of  the  nominal  value. 


5.3  Conclusion 

A  variational  formulation  for  the  depth  problem  of  a  shallow  water  Pekeris  waveguide  has  been 
presented  and  eigensolution  obtained  using  a  Rayleigh-Ritz  approximation,  for  low  to  medium 
frequency  range.  An  interesting  feature  of  the  model  is  that  the  trial  functions  are  derived  from  an 
isovelocity  problem  that  has  exact  solution.  It  is  important  to  note  that  such  trial  functions 
automatically  satisfy  even  the  dynamic  interface  condition  at  the  seabed,  thus  contributing  to  the 
accuracy  of  the  numerical  model.  The  model  has  been  implemented  in  MATLAB.  The  radial 
wavenumbers  obtained  compare  very  well  with  KRAKEN  results.  The  proposed  model  is  accurate 
and  is  deemed  to  provide  a  compact  fonn  for  the  depth  modes  as  a  linear  combination  of  the 
isovelocity  modes. 

Admittedly,  the  computational  effort  in  setting  up  the  matrix  in  the  proposed  RR  model  using 
numerical  quadrature  is  high  compared  to  setting  up  the  finite  difference  based  matrix  in  the 
KRAKEN  code.  However,  noting  the  diagonal  dominance  of  the  matrix  obtained  in  the  RR  model, 
it  would  be  worthwhile  exploring  the  possibility  of  approximating  it  by  a  narrow  banded  matrix  in 
order  to  reduce  the  volume  of  computation  in  setting  up  the  matrix  and  possibly  in  obtaining  the 
eigensolution. 

The  proposed  RR  model  presented  readily  admits  heterogeneous  fluid  halfspace,  although  only 
isovelocity  cases  have  been  studied.  Extension  of  the  Rayleigh-Ritz  model  to  high  frequency 
shallow  water  waveguides  as  well  as  medium  and  high  frequency  deep  water  waveguides  would 
seem  possible  in  principle.  However,  since  such  problems  will  have  a  few  to  several  hundred 
modes,  the  perfonnance  of  the  RR  model  remains  to  be  proved.  Extension  to  shallow  water 
geoacoustic  waveguides  would  be  worth  investigating. 
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Figure  5.2(a)  Sound  speed  profile- 1  Figure  5.2(b)  Sound  speed  profile-2  (Schneider,  1977) 
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Figure  5.2c  Sound  speed  profde-3  (Sazontov  et  al,  2002) 
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Chapter  6 
Epilogue 


6.1  Summary 

The  main  thrust  of  this  study  has  been  to  fonnulate  a  finite  element  model  for  elastic  scattering  by 
objects  in  a  range  and  depth  dependent  cylindrical  ocean  acoustic  waveguide.  The  3-D  wave 
equation  for  acoustic  pressure  has  been  chosen.  Assuming  time-harmonic  motion,  the  acoustic 
field  in  the  azimuthal  direction  has  been  expressed  in  a  Fourier  series,  leading  to  a  2-D  Helmholtz 
equation  for  each  of  the  hannonics.  The  acoustic  field  due  to  linear  elastic  scattering  by  an 
axisymmetric  body  in  the  waveguide  has  been  expressed  as  a  sum  of  rigid-scattered  and  radiated 
acoustic  fields.  Hannonic  acoustic  radiation  by  an  axisymmetric  body  has  been  discussed 
employing  an  assumed  mode  approach  to  structural  displacements.  The  acoustic  radiation 
(pressure)  loading  on  the  structure  has  been  identified  in  terms  of  added  mass  and  radiation 
damping  matrices,  which  can  be  appended  to  the  structural  mass  and  damping  matrices,  thus 
completing  the  fluid-structure  interaction  formulation.  The  acoustic  loading  due  to  rigid-scattered 
field  constitutes  the  forcing  on  the  coupled  system. 

A  finite  element  model  has  been  fonnulated  to  solve  the  2-D  Helmholtz  equation  for  each  of  the 
circumferential  harmonics  by  truncating  the  fluid  domain  at  a  finite  range.  A  Bayliss-type 
boundary  operator  has  been  adopted  for  each  of  the  waveguide  propagating  modes,  to  serve  as  an 
approximate  nonreflecting/absorbing  condition  on  the  truncation  boundary.  The  proposed  FE 
model  applies  to  the  empty  waveguide  with  a  point  source,  the  waveguide  embedded  with  a  rigid 
scattering  object  and  the  waveguide  with  the  scattering  object  considered  as  a  radiator.  The 
generalized  FE  model  has  been  implemented  in  a  FORTRAN90  code. 

Depth  eigensolution  of  the  depth  dependent  waveguide  is  required  to  impose  the  nonreflecting 
boundary  condition  at  the  truncation  boundary  in  the  FE  model.  A  Rayleigh-Ritz  approximation 
has  been  developed  for  this  purpose,  using  the  depth  eigenmodes  of  the  isovelocity  problem.  The 
RR-model  has  been  implemented  using  MATLAB  coding. 

6.2  Concluding  Remarks 

The  FSI  formulation  for  harmonic  elastic  scattering  by  axisymmetric  objects  in  a  waveguide 
conveniently  decouples  the  acoustic/fluid  dynamic  formulation  from  the  structural  formulation, 
the  former  being  the  primary  focus  of  this  study.  The  FE  model  for  the  2-D  Helmholtz  equation, 
despite  being  computationally  intensive,  readily  accounts  for  complex  axisymmetric  geometries, 
multiple  scattering  and  effect  of  free  surface  and  seafloor. 

It  is  interesting  to  note  in  the  present  study  while  an  elaborate  finite  element  model  has  been 
adopted  for  representing  heterogeneous  waveguide  propagation,  a  classical  approach  has  been 
indicated  for  modeling  the  vibration  of  axisymmetric  scatterers,  instead  of  a  finite  element  model. 
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Despite  its  limitations,  it  is  deemed  that  the  structural  fonnulation  suggested  provides  a  compact 
model,  which  in  turn  simplifies  the  finite  element  model  for  the  waveguide  radiation  problem. 

The  Rayleigh-Ritz  approximation  developed  for  the  depth  eigenproblem  yields  eigenmodes  as  a 
linear  combination  of  isovelocity  modes.  Such  a  form  is  more  compact  than  the  eigenmodes 
derived  from  finite  difference  or  finite  element  models  for  the  depth  problem,  and  hence 
convenient  to  use  with  the  FE  model  for  the  waveguide. 

6.3  Further  Research 

The  FORTRAN  code  developed  in  this  study  is  fairly  large  in  size  and  extensive  testing  has  been 
done  to  validate  different  components  of  the  numerical  model.  FE  meshing  and  data  generation 
code  has  been  developed  to  exploit  a  commercially  available  standard  FE  preprocessor.  As  of 
now,  coding  for  the  radiation  problem  is  complete  and  validation  examples  remain  to  be  solved. 
The  scattering  problem  involves  extracting  nonaxisymmetric  scattered  pressure  on  the  surface  of 
an  axisymmetric  body.  Coding  for  this  task  is  in  progress.  The  research  work  reported  here  fonns 
part  of  the  doctoral  work  of  the  second  author  of  this  report.  So,  all  the  remaining  tasks  will  be 
brought  to  a  logical  conclusion  and  included  in  the  doctoral  thesis. 
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Table  5.1  Propagating  mode  wavenumbers  for  100  m  deep  single  layer  waveguide 
with  pressure  release  top  and  rigid  bottom,  for  sound  speed  profile- 1  in  Fig.  5.2a 


Frequency  50  Hz 

Frequency  100  Hz 

Frequency  500  Hz 

Frequency  1000  Hz 

j 

Kraken 

RR  Model 

(n  =  2  np) 

j 

Kraken 

RR  Model 

(n  =  2  np) 

j 

Kraken 

RR  Model 

(n  =  2  np) 

j 

Kraken 

RR  Model 

(«  =  2  np) 

1 

0.20626 

214 

0.20626 

201 

1 

0.41370 

960 

0.41370 

970 

1 

2.07327 

564 

2.07327 

560 

1 

4.14861 

035 

4.14861 

026 

2 

0.20034 

991 

0.20034 

985 

2 

0.40927 

711 

0.40927 

704 

2 

2.06845 

466 

2.06845 

461 

2 

4.14127 

043 

4.14127 

031 

3 

0.19017 

967 

0.19017 

943 

3 

0.40377 

782 

0.40377 

775 

3 

2.06486 

547 

2.06486 

537 

3 

4.13829 

130 

4.13829 

125 

4 

0.17395 

532 

0.17395 

521 

6 

0.37362 

274 

0.37362 

269 

29 

1.85314 

779 

1.85314 

761 

64 

3.60056 

038 

3.60056 

Oil 

5 

0.14963 

467 

0.14963 

417 

7 

0.35742 

699 

0.35742 

694 

30 

1.83761 

895 

1.83761 

876 

65 

3.58298 

492 

3.58298 

469 

6 

0.11185 

534 

0.11185 

467 

8 

0.33754 

919 

0.33754 

916 

31 

1.82146 

158 

1.82146 

144 

66 

3.56505 

546 

3.56505 

524 

7 

0.02586 

682 

0.02586 

071 

11 

0.24621 

631 

0.24621 

595 

64 

0.50679 

833 

0.50679 

799 

129 

0.80570 

636 

0.80570 

489 

12 

0.19734 

847 

0.19734 

731 

65 

0.36129 

618 

0.36129 

551 

130 

0.62813 

151 

0.62812 

972 

13 

0.12350 

950 

0.12350 

583 

66 

0.04717 

033 

0.04716 

571 

131 

0.37142 

513 

0.37142 

194 

Table  5.2  Propagating  mode  wavenumbers  for  100  m  deep  single  layer  waveguide 
with  pressure  release  top  and  rigid  bottom,  for  sound  speed  profile-2  in  Fig.  5.2b  [10] 


Frequency  50  Hz 

Frequency  100  Hz 

Frequency  500  Hz 

Frequency  1000  Hz 

j 

Kraken 

RR  Model 

( n  =  2  np) 

j 

Kraken 

RR  Model 

(n  =  2  np) 

j 

Kraken 

RR  Model 

(n  =  2  np) 

j 

Kraken 

RR  Model 

(n  =  2  np) 

1 

0.21750 

372 

0.21750 

357 

1 

0.43587 

227 

0.43587 

231 

1 

2.18221 

704 

2.18221 

738 

1 

4.36602 

483 

4.36602 

623 

2 

0.21290 

853 

0.21290 

720 

2 

0.43359 

244 

0.43359 

229 

2 

2.18053 

321 

2.18053 

347 

2 

4.36375 

413 

4.36375 

598 

3 

0.20324 

210 

0.20323 

918 

3 

0.42871 

517 

0.42871 

464 

3 

2.17934 

145 

2.17934 

179 

3 

4.36201 

332 

4.36201 

389 

4 

0.18792 

005 

0.18791 

626 

6 

0.39951 

802 

0.39951 

882 

34 

1.90610 

439 

1.90610 

477 

69 

3.78584 

771 

3.78584 

678 

5 

0.16548 

007 

0.16547 

682 

7 

0.38456 

842 

0.38456 

818 

35 

1.88843 

941 

1.88843 

974 

70 

3.76780 

058 

3.76779 

971 

6 

0.13237 

311 

0.13237 

254 

8 

0.36632 

398 

0.36632 

402 

36 

1.87008 

454 

1.87008 

450 

71 

3.74945 

258 

3.74945 

203 

7 

0.07556 

731 

0.07556 

539 

12 

0.24303 

274 

0.24302 

958 

67 

0.61368 

851 

0.61368 

978 

137 

0.75886 

488 

0.75885 

584 

13 

0.18816 

133 

0.18815 

815 

68 

0.49430 

767 

0.49430 

752 

138 

0.55264 

353 

0.55263 

741 

14 

0.09895 

384 

0.09894 

815 

69 

0.33193 

990 

0.33194 

354 

139 

0.18167 

275 

0.18165 

878 
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Table  5.3  Propagating  mode  wavenumbers  for  150  m  deep  single  layer  waveguide 
with  pressure  release  top  and  rigid  bottom,  for  sound  speed  profile-3  [11] 


Frequency  50  Hz 

Frequency  100  Hz 

Frequency  500  Hz 

Frequency  1000  Hz 

j 

Kraken 

RR  Model 

0 n  =  2  nP) 

j 

Kraken 

RR  Model 

(n  =  2  np) 

j 

Kraken 

RR  Model 
(n  =  2  rip) 

j 

Kraken 

RR  Model 

(n  =  2np) 

1 

0.21379 

035 

0.21379 

031 

1 

0.42807 

733 

0.42807 

730 

1 

2.14182 

484 

2.14182 

465 

1 

4.28393 

384 

4.28393 

346 

2 

0.21118 

770 

0.21118 

766 

2 

0.42618 

778 

0.42618 

772 

2 

2.14065 

208 

2.14065 

178 

2 

4.28324 

895 

4.28324 

844 

3 

0.20679 

341 

0.20679 

334 

3 

0.42346 

331 

0.42346 

321 

3 

2.13887 

874 

2.13887 

849 

3 

4.28158 

839 

4.28158 

770 

4 

0.20045 

527 

0.20045 

521 

8 

0.39690 

255 

0.39690 

250 

49 

1.87715 

622 

1.87715 

607 

100 

3.72551 

546 

3.72551 

528 

5 

0.19153 

272 

0.19153 

267 

9 

0.38798 

638 

0.38798 

633 

50 

1.86567 

414 

1.86567 

401 

101 

3.71372 

757 

3.71372 

738 

6 

0.17966 

348 

0.17966 

341 

10 

0.37766 

117 

0.37766 

108 

51 

1.85388 

898 

1.85388 

886 

102 

3.70177 

828 

3.70177 

810 

7 

0.16438 

144 

0.16438 

129 

11 

0.36586 

465 

0.36586 

456 

52 

1.84178 

203 

1.84178 

192 

103 

3.68967 

655 

3.68967 

636 

8 

0.14451 

594 

0.14451 

157 

18 

0.21885 

397 

0.21885 

275 

100 

0.46157 

064 

0.46157 

029 

202 

0.64265 

825 

0.64265 

713 

9 

0.11775 

220 

0.11775 

171 

19 

0.17917 

147 

0.17916 

949 

101 

0.35400 

344 

0.35400 

296 

203 

0.48558 

547 

0.48558 

401 

10 

0.07725 

818 

0.07725 

680 

20 

0.12424 

114 

0.12423 

713 

102 

0.19161 

962 

0.19161 

871 

204 

0.24022 

424 

0.24022 

120 

Table  5.4  Propagating  mode  wavenumbers  for  two  layer  waveguide  -  top  layer  100  m  deep  with  sound  speed  profile- 1;  isovelocity 

bottom  layer  50  m  deep(  p  =  1024  kg/  m3 ;  p2  =  2000  0  kg/  m3 ;  c2  =  2000  m/s) 


Frequency  50  Hz 

Frequency  100  Hz 

Frequency  500  Hz 

Frequency  1000  Hz 

j 

Kraken 

RR  Model 

0 n  =  2  np) 

j 

Kraken 

RR  Model 

(n  =  2  np) 

j 

Kraken 

RR  Model 

(«  =  2  np) 

j 

Kraken 

RR  Model 

{n  =  2  np) 

1 

0.20456 

112 

0.20456 

103 

1 

0.41203 

261 

0.41203 

272 

1 

2.07030 

231 

2.07030 

222 

1 

4.14439 

339 

4.14439 

290 

2 

0.19822 

067 

0.19822 

067 

2 

0.40748 

072 

0.40748 

071 

2 

2.06729 

366 

2.06729 

348 

2 

4.13928 

774 

4.13928 

750 

3 

0.18785 

719 

0.18785 

687 

3 

0.40192 

320 

0.40192 

319 

3 

2.06292 

091 

2.06292 

079 

3 

4.13716 

472 

4.13716 

405 

4 

0.17204 

366 

0.17204 

356 

7 

0.35523 

736 

0.35523 

732 

45 

1.56626 

702 

1.56627 

064 

90 

3.12889 

039 

3.12889 

046 

5 

0.15456 

298 

0.15456 

313 

8 

0.33584 

249 

0.33584 

244 

46 

1.56266 

492 

1.56266 

466 

91 

3.12285 

967 

3.12285 

252 

B 

0.14850 

700 

0.14850 

645 

9 

0.31364 

127 

0.31364 

508 

47 

1.55526 

091 

1.55526 

088 

92 

3.11923 

353 

3.11924 

959 

7 

0.12763 

872 

0.12763 

890 

16 

0.19557 

802 

0.19557 

645 

89 

0.38141 

305 

0.38141 

320 

179 

0.62752 

845 

0.62753 

139 

8 

0.10837 

090 

0.10836 

998 

17 

0.14851 

788 

0.14851 

899 

90 

0.29293 

460 

0.29293 

311 

180 

0.47294 

187 

0.47293 

210 

9 

0.05543 

307 

0.05543 

244 

18 

0.10714 

761 

0.10713 

932 

91 

0.04212 

166 

0.04212 

181 

181 

0.32370 

017 

0.32370 

695 
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Table  5.5  Wavenumbers  for  Pekeris  waveguide  example  at  100  Hz  [9] 

( p  =  1000  kg/  m 3 ;  c  =  1500  m/s  ;  ph  =  p2  =  1250  kg/  m 3 ;  c2  =  1600  m/s; Attenuation  in  half  space  =  0.3 125  db/km-Hz) 


j 

Kraken 

Ref.  9 

RR  Model(n  =  2np) 

1 

0.41585  560 

+  i  0.54834E-06 

0.41585  35 

+  iO. 0000547 

0.41585  345 

+  i0.54674E-06 

2 

0.40668  711 

+  i  0.24050E-05 

0.40668  12 

+  iO. 0002392 

0.40668  125 

+  i  0.23923E-05 

3 

0.39183  488 

+  iO.  13149E-04 

0.39184  12 

+  iO. 0013215 

0.39184  127 

+  i  0.13216E-04 
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Table  5.6  Propagating  mode  wave  numbers  for  100  m  deep  PEKERIS  waveguide-  sound  speed  profile- 1  in  water  column 


( p  =  1024  kg/ m3 ;  ph  =  p2  =2000  0  kg/ m3 ;  c2  =  2000  m/s;  Attenuation  =0.3 125  db/km-Hz )  [9] 


Frequency  50  Hz 

Frequency  100  Hz 

Frequency  500  Hz 

Frequency  1000  Hz 

j 

Kraken 

RR  Model 

0  =  2  np) 

j 

Kraken 

RR  Model 

(n  =  2  np) 

j 

Kraken 

RR  Model 

( n  =  2  np) 

j 

Kraken 

RR  Model 

(n  =  2np) 

0.20456  091 

0.20456  096 

0.41203  246 

0.41203  255 

2.07030  219 

2.07030  211 

4.14439  324 

4.14439  277 

1 

+ 

+ 

1 

+ 

+ 

1 

+  i 

+  i 

1 

+  i 

+  i 

i0.10143E-06 

iO.101 1  IE-06 

i0.63370E-07 

i0.63402E-07 

0.4951  IE-07 

0.4545E-07 

0.61407E-07 

0.61452E-07 

0.19822  000 

0.19822  053 

0.40748  043 

0.40748  040 

2.06729  356 

2.06729  339 

4.13928  765 

4.13928  742 

2 

+ 

i0.25266E-06 

+  i0.2536E-06 

2 

+ 

i0.12214E-06 

+ 

i0.12227E-06 

2 

+ 

i0.40969E-07 

+  i 

0.40993E-07 

2 

+  i 

0.35846E-07 

+  i 

0.35863E-07 

0.18785  560 

0.18785  605 

0.40192  273 

0.40192  272 

2.06292  074 

2.06292  064 

4.13716  457 

4.13716  395 

3 

+ 

i0.48573E-06 

+  i0.4901E-06 

3 

+ 

iO.  1 861 8E — 06 

+ 

i0.18642E-06 

3 

+  i 

0.6733  IE-07 

+  i 

0.67385E-07 

3 

+  i 

0.56470E-07 

+  i 

0.56524E-07 

0.17203  622 

0.17204  994 

0.39435  248 

0.39435  260 

1.96218  955 

1.96218  975 

3.91178  541 

3.91178  556 

4 

+ 

+  i0.9934E-06 

4 

+ 

+ 

20 

+  i 

+  i 

41 

+  i 

+  i 

i0.96521E-06 

i0.27286E-06 

i0.27324E-06 

0.31449E-06 

0.31528E-06 

0.33135E-07 

0.33219E-06 

0.38427  786 

0.38427  807 

1.95219  366 

1.95219  396 

3.90148  034 

3.90148  040 

5 

+ 

+ 

21 

+  i 

+  i 

42 

+  i 

+  i 

i0.40048E-06 

i0.401 10E-06 

0.33974E-06 

0.34062E-06 

0.34428E-06 

0.34519E-06 

0.37125  268 

0.37125  300 

1.94159  584 

1.94159  617 

3.89085  762 

3.89085  781 

6 

+ 

+ 

22 

+  i 

+  i 

43 

+  i 

+  i 

i0.56060E-06 

i0.56234E-06 

0.36500E-06 

0.36587E-06 

0.35895E-06 

0.35989E-06 

0.35523  437 

0.35523  489 

1.64060  761 

1.64061  287 

3.19707  658 

3.19709  962 

7 

+ 

+ 

40 

+  i 

+  i 

83 

+  i 

+  i 

i0.79275E-06 

i0.79506E-06 

0.17135E-05 

0.17410E-05 

0.28194E-06 

0.28947E-05 

0.33583  530 

0.33583  538 

1.61662  367 

1.61663  562 

3.17155  536 

3.17159  690 

8 

+ 

+ 

41 

+  i 

+  i 

84 

+ 

+  i 

i0.12559E-05 

i0.12388E-05 

0.21630E-05 

0.21800E-05 

i0.37810E-05 

0.38429E-05 

1.59175  377 

1.59177  107 

3.14553  001 

3.14554  474 

42 

+  i 

+  i 

85 

+  i 

+  i 

0.32057E-05 

0.32433E-05 

0.77887E-05 

0.76853E-05 
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Appendix 


The  Rayleigh-Ritz  model  presented  for  the  depth  eigenproblems,  discussed  in  Sec.  5.3,  employs 
the  analytical  depth  modes  for  an  isovelocity  waveguide  as  the  trial  functions.  The  details  of  the 
various  isovelocity  waveguide  examples  used  are  presented  here. 


For  a  single  layer  waveguide  shown  in  Fig.  A.  1,  the  trial  functions  are  given  by  [1,  2,  5] 

,  (/-0.5W 

Wj  =  sin (kZJz)  ’  where  k--i  =  - =  1’2’- 

For  a  two  layer  waveguide  shown  in  Fig.  A.2,  the  trial  functions  are  given  by  [1] 
y/i  (z)  =  a  sin  kzUz  Q<z<Dx 


where  kzXi 


_a  sin^i  A 
cos  £z2i(A-A) 
( o1lcl-k;,k:2i=co1l 


cos  kz2i{D2-z)  Dx<  z<D2 

cl  -  kl  and  kr  is  the  solution  of  the  transcendental  equation 


k  o 

2,1  2  =  tan  k  A  tan  k  0h1 

i  Zli  1  Zl  A  Z 

KiPi 


For  a  Pekeris  waveguide  (for  which  D2  -^  qo  in  Fig.  A.2  the  trial  functions  are  given  by  [1,  2] 
y/t  (z)  =  a  sin  kzUz  0  <  z  <  D{ 

=  a  sin  kzlih  e~‘kz2,<-h~z)  Dx  <  z  <  oo 


Figure  A.  1  Single  layer  waveguide 


Figure  A.2  Two  layer  waveguide 
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